
files = ['RL';'RR';'LL';'LR'];
Ts = [14,20,30,40,50,55,60,65,70,80,100,130,160,190,250,296];
Rcoeffs=[1.133,1.121];
offsets=[]
%Ts = [80]%,100,130,160,190,220,250,275,296];
figure
hold on
EsVH=zeros(16,991);
dataVH=zeros(16,991);
for num=1:length(Ts)
    filename=['LL',num2str(Ts(num)),'K.txt'];
        data=readmatrix(filename);
        plot(data(26:end,1),data(26:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        dataVH(num,:)=data(26:(26+990),2);
        EsVH(num,:)=data(26:(26+990),1);
        ylim([20,1000])
        title('LL')
end   

figure
[X,Y] = meshgrid(EsVH(1,:),Ts);
pcolor(X,Y,dataVH)

shading interp
colormap('jet')
clim([10,50])

pol = ['RL';'RR';'LL';'LR'];
%Ts = [80]%,100,130,160,190,220,250,275,296];
figure
hold on
EsHV=zeros(16,991);
dataHV=zeros(16,991);
for num=1:length(Ts)
    filename=['RR',num2str(Ts(num)),'K.txt'];
        data=readmatrix(filename);
        plot(data(26:end,1),data(26:end,2)+50*(num-1),'.-','Color',[(num)/20 0 0]);
        dataHV(num,:)=data(26:(26+990),2);
        EsHV(num,:)=data(26:(26+990),1);
        ylim([20,1000])
        title('RR')
end
figure
[X,Y] = meshgrid(EsHV(1,:),Ts);
    pcolor(X,Y,dataHV)

shading interp
colormap('jet')
clim([10,200])

%% offset2

offsetsVH=zeros(16,1);
offsetsHV=zeros(16,1);
offseteqn='a1*exp(-(x-b1)^2/c1^2)+d1';
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',[10001,0,2,10]);
ft = fittype(offseteqn,'options',fo);


for num=1:length(Ts)
    EsVHfit2=EsVH(num,:);
    dataVHfit2=dataVH(num,:);
    EsHVfit2=EsHV(num,:);
    dataHVfit2=dataHV(num,:);
    EsHVfit2(dataHVfit2>10000)=300;
    dataHVfit2(dataHVfit2>10000)=0;
    EsVHfit2(dataVHfit2>10000)=300;
    dataVHfit2(dataVHfit2>10000)=0;
    [curve2,gof2] = fit(EsVHfit2',dataVHfit2',ft);
    offsetsVH(num)=curve2.b1;
    figure(22)
    plot(EsVHfit2',dataVHfit2','o')
    
    hold on
    plot(-10:0.01:10,curve2(-10:0.01:10))
    xlim([-10,10])
    hold off
    figure(23)
    [curve3,gof2] = fit(EsHVfit2',dataHVfit2',ft);
    plot(EsHVfit2',dataHVfit2','o')
    
    hold on
    plot(-10:0.01:10,curve3(-10:0.01:10))
    xlim([-10,10])
    hold off
    offsetsHV(num)=curve3.b1;
end

offsets= offsetsHV/2+offsetsVH/2;


%% Gaussfit

EsHVfit=EsHV(:,700:800);
EsVHfit=EsVH(:,700:800);

dataHVfit=dataHV(:,700:800);
dataVHfit=dataVH(:,700:800);
dataHVfit(11,35)=dataVHfit(11,35);
dataHVfit(11,36)=dataVHfit(11,36);
eqn='a1*exp(-(x-b1)^2/c1^2)+d1';
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',[20,298,10,10]);
ft = fittype(eqn,'options',fo);
asHV=[]
bsHV=[]
bsHVerr=[]
dsHV=[]
figure
hold on
for idx= 1:length(Ts)
    E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    asHV=[asHV,curve2.a1];
    bsHV=[bsHV,curve2.b1];
    dsHV=[dsHV,curve2.d1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr,berr(2,2)/2-berr(1,2)/2];
    plot(E,I+30*(idx-1),'.-','Color',[Ts(idx)/300 0 0])
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1),'Color',[Ts(idx)/300 0 0],'LineWidth',1.5)
end
xlim([285,310])
title('RR')
figure(30);
errorbar(Ts,bsHV-mean(offsetsHV'),bsHVerr,'o')
asVH=[]
bsVH=[]
bsVHerr=[]
dsVH=[]
figure
hold on
for idx= 1:length(Ts)
    E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    asVH=[asVH,curve2.a1];
    bsVH=[bsVH,curve2.b1];
    dsVH=[dsVH,curve2.d1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr,berr(2,2)/2-berr(1,2)/2];
    plot(E,I+30*(idx-1),'.-','Color',[Ts(idx)/300 0 0])
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1),'Color',[Ts(idx)/300 0 0],'LineWidth',1.5)
end
xlim([285,310])
title('LL')
figure(30);
hold on
errorbar(Ts,bsVH-mean(offsetsVH'),bsVHerr,'o')

figure();
hold on
plot(Ts,asHV,'o')
plot(Ts,asVH,'o')
%% imagesc plot
figure
[X,Y] = meshgrid([EsHV(1,700:800)],[Ts]);

pcolor(X-mean(offsets),Y,[dataHVfit-dsHV'+dataVHfit-dsVH']/2)

shading interp
colormap jet
clim([0,40])
xlim([285,310])
ylim([0,300])
title('LL+RR532')


figure
[X,Y] = meshgrid(EsVH(1,:),Ts);
pcolor(X-mean(offsetsVH),Y,dataVH-dsVH')
shading interp
colormap jet
clim([0,40])
xlim([285,310])
ylim([0,300])


Tplot=0:305;
dataplot=zeros(length(Tplot),991*2);
offsetplotHV=[];
offsetplotVH=[];
Tline=[];
for num=1:(length(Ts)-1)
    Tline(num)=(Ts(num)/2+Ts(num+1)/2)-0.1;
end
for idx=1:length(Tplot)
    numT=sum(Tplot(idx)>Tline)+1
    dataplot(idx,:)=[dataHV(numT,:),dataVH(numT,:)];
    offsetplotHV=[offsetplotHV,offsetsHV(numT)];
    offsetplotVH=[offsetplotVH,offsetsVH(numT)];
end


figure
[X,Y] = meshgrid([EsHV(1,:),EsVH(1,:)],Tplot);
pcolor([X(:,1:991)-offsetplotHV',X(:,992:1982)-offsetplotVH'],Y,dataplot)

shading interp
colormap jet
xlim([285,310])
clim([10,50])
%% 
figure(30)
hold on
bsall1=[bsVH-mean(offsetsVH'),bsHV-mean(offsetsHV')];
Tsall1=[Ts,Ts];
a=1.834;
c=299.7;
a=1.9174;
c=299.8391;
eqnfit=@(x) -a*(1+2./(exp((c)./(2*0.69501*x))-1))+c;
plot(0:300,eqnfit(0:300))
a=1.81;
c=299.7;
a=2.0069;
c=300.0777;
eqnfit=@(x) -a*(1+2./(exp((c)./(2*0.69501*x))-1))+c;
bspart1=[bsVH(9:end)-mean(offsetsVH(9:end)'),bsHV(9:end)-mean(offsetsHV(9:end)')];
Tspart1=[Ts(9:end),Ts(9:end)];
plot(0:300,eqnfit(0:300))
%ylim([290,300])

%eqnfit1='-a*(1+2/(exp((c)/(2*0.69501*x))-1))+c';

%%
%{
bsall=[bsall1,bsall2];
Tsall=[Tsall1,Tsall2];
bspart=[bspart1,bspart2];
Tspart=[Tspart1,Tspart2];
a=1.9174;
c=299.8391;
a=2.0069;
c=300.0777;
%}
%% seperate peaks
bs=[293,295,297,299,300.9,303];
eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+e*x+f';
sp=[10,30,50,30,10, bs(1), bs(2), bs(3), bs(4), bs(5),.4,0,0];
eqn2='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+e*x+f';
sp2=[10,10,10, bs(2)-1.1, bs(3)-1.1, bs(4)-1.1,1,0,0];
fo = fitoptions('Method','NonlinearLeastSquares',...
            'StartPoint',sp);
ft = fittype(eqn,'options',fo);
bsHV=[]
bsHVerr=[]
dsHV=[]
figure
hold on
for idx= 1:length(Ts)
    if Ts(idx)<68
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    
    E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsHV=[bsHV;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5]];
    dsHV=[dsHV,curve2.c1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr;[berr(2,6)/2-berr(1,6)/2, berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    else
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp2);
        ft = fittype(eqn2,'options',fo);
        E=EsHVfit(idx,:);
    I=dataHVfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsHV=[bsHV;[0,curve2.b1,curve2.b2,curve2.b3,0]];
    dsHV=[dsHV,curve2.c1];
    berr=confint(curve2,0.68);
    bsHVerr=[bsHVerr;[0, berr(2,4)/2-berr(1,4)/2, berr(2,5)/2-berr(1,5)/2, berr(2,6)/2-berr(1,6)/2,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
end

figure(40);
errorbar(Ts,bsHV-mean(offsetsHV),bsHVerr,'o','Color',[0.85,0.33,0.10],'LineWidth',1)
ylim([290,310])

bsVH=[]
bsVHerr=[]
dsVH=[]
figure
hold on
for idx= 1:length(Ts)
    if Ts(idx)<68
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    
    E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsVH=[bsVH;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5]];
    dsVH=[dsVH,curve2.c1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr;[berr(2,6)/2-berr(1,6)/2, berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    else
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp2);
        ft = fittype(eqn2,'options',fo);
        E=EsVHfit(idx,:);
    I=dataVHfit(idx,:);
    [curve2,gof2] = fit(E',I',ft);
    bsVH=[bsVH;[0,curve2.b1,curve2.b2,curve2.b3,0]];
    dsVH=[dsVH,curve2.c1];
    berr=confint(curve2,0.68);
    bsVHerr=[bsVHerr;[0, berr(2,4)/2-berr(1,4)/2, berr(2,5)/2-berr(1,5)/2, berr(2,6)/2-berr(1,6)/2,0]];
    plot(E,I+30*(idx-1),'.-')
    plot(285:0.1:310, curve2(285:0.1:310)+30*(idx-1))
    end
end
figure(43);
hold on
errorbar(Ts,bsVH-mean(offsetsVH),bsVHerr,'diamond','Color',[0.93,0.69,0.13],'LineWidth',1)

ylim([290,306])
xlim([0,110])


figure(41);
hold on
plot(Ts,dsVH*2,'diamond','Color',[0.93,0.69,0.13],'LineWidth',1)
plot(Ts,dsHV*2,'o','Color',[0.85,0.33,0.10],'LineWidth',1)


ts= [Ts(1:11),Ts(1:11)]
ds= [dsVH(1:11),dsHV(1:11)]



%% seperate peaks LL+RR
bs=[293,295,297,299,300.9,303];
eqn='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+a4*exp(-(x-b4)^2/c1^2)+a5*exp(-(x-b5)^2/c1^2)+e*x+f';
sp=[19,61,75,65,30, bs(1)-.2, bs(2)-.2, bs(3)-.2, bs(4)-.2, bs(5)-.2,.4,0.02,28];
eqn2='a1*exp(-(x-b1)^2/c1^2)+a2*exp(-(x-b2)^2/c1^2)+a3*exp(-(x-b3)^2/c1^2)+e*x+f';
sp2=[10,10,10, bs(2)-2, bs(3)-2, bs(4)-2,0.4,0,0];
eqn3='a1*exp(-(x-b1)^2/c1^2)+e*x+f';
sp3=[10, bs(3)-2,5,0,0];

bs = []; 
bserr = [];
bsHV=[]
bsHVerr=[]
dsHV=[]
dsHVerr=[]
figure
hold on
for idx= 1:length(Ts)
    fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp);
        ft = fittype(eqn,'options',fo);
    if Ts(idx)<102
        E=EsHVfit(idx,6:66);
        I=dataHVfit(idx,6:66)+dataVHfit(idx,6:66);
        [curve2,gof2] = fit(E',I',ft);
        bsHV=[bsHV;[curve2.b1,curve2.b2,curve2.b3,curve2.b4,curve2.b5]];
        dsHV=[dsHV,curve2.c1];
        berr=confint(curve2,0.68);
        bsHVerr=[bsHVerr;[berr(2,6)/2-berr(1,6)/2, berr(2,7)/2-berr(1,7)/2, berr(2,8)/2-berr(1,8)/2, berr(2,9)/2-berr(1,9)/2, berr(2,10)/2-berr(1,10)/2]];
        dsHVerr=[dsHVerr;[berr(2,11)/2-berr(1,11)/2]];
        plot(E,I+80*(idx-1),'.','Color',[min([Ts(idx)/100,1]) 0 0])
        plot(285:0.1:310, curve2(285:0.1:310)+80*(idx-1),'Color',[min([Ts(idx)/100,1]) 0 0])
    end

    if Ts(idx)>102
        fo = fitoptions('Method','NonlinearLeastSquares',...
                'StartPoint',sp3);
        ft = fittype(eqn3,'options',fo);
        E=EsHVfit(idx,6:66);
        I=dataHVfit(idx,6:66)+dataVHfit(idx,6:66);
        [curve2,gof2] = fit(E',I',ft);
        %bsHV=[bsHV;[0,curve2.b1,curve2.b2,curve2.b3,0]];
        %dsHV=[dsHV,curve2.c1];
        berr=confint(curve2,0.68);
        %bsHVerr=[bsHVerr;[0, berr(2,4)/2-berr(1,4)/2, berr(2,5)/2-berr(1,5)/2, berr(2,6)/2-berr(1,6)/2,0]];
        %dsHVerr=[dsHVerr;[berr(2,7)/2-berr(1,7)/2]];
        plot(E,I+80*(idx-1),'.','Color',[min([Ts(idx)/100,1]) 0 0])
        plot(285:0.1:310, curve2(285:0.1:310)+80*(idx-1),'Color',[min([Ts(idx)/100,1]) 0 0])
    end
    I=I-E*curve2.e-curve2.f;
    weighted_avg = sum(E .* I) / sum(I);
    
    % Compute standard error
    weighted_var = sum((E - weighted_avg).^2 .* I) / sum(I);
    weighted_err = sqrt(weighted_var / length(I));
    
    % Store results
    bs = [bs, weighted_avg];
    bserr = [bserr, weighted_err];
end



figure(42);
hold on
errorbar(Ts,bs,bserr,'o','Color','r','LineWidth',1,'MarkerSize',4)
ylim([290,304])



%% Weighted sum LL+RR

bs = []; 
bserr = [];

for idx = 1:length(Ts)
    fo = fitoptions('Method', 'NonlinearLeastSquares', ...
                    'StartPoint', sp);
    ft = fittype(eqn, 'options', fo);
    
    E = EsHVfit(idx, :);
    I = dataHVfit(idx, :) + dataVHfit(idx, :);
    
    % Compute weighted average
    weighted_avg = sum(E .* I) / sum(I);
    
    % Compute standard error
    weighted_var = sum((E - weighted_avg).^2 .* I) / sum(I);
    weighted_err = sqrt(weighted_var / length(I));
    
    % Store results
    bs = [bs, weighted_avg];
    bserr = [bserr, weighted_err];
end

figure(43);
hold on
errorbar(Ts,bs,bserr,'o','Color','r','LineWidth',1,'MarkerSize',4)
ylim([290,304])


figure(41);
hold on
errorbar(Ts,dsVH,dsHVerr,'o','Color','r','LineWidth',1,'MarkerSize',4)
xlim([0,110])
ylim([0.4,1])





















